% Plot model theta against pi using the two-productivity model
% October 2022
% Makoto Nirei

clearvars;
setparams;

pilength=100;
piv=linspace(0.0005,0.5,pilength);
theta_v=zeros(1,pilength); w_v=zeros(1,pilength); pb_v=zeros(2,pilength);
for ind_pi=1:pilength
    pi=piv(ind_pi);
    
    get_equilibrium;
    
    theta=0;
    for i=1:2
        theta = theta + 0.5* pb(i)^(1-eta) * vp(q(i),1-eta)/vp(q(i),mu/pi);
    end
    theta_v(ind_pi)=theta;
    w_v(ind_pi)=w;
    pb_v(:,ind_pi)=reshape(pb,2,1);
end

figure(3),
plot(piv,theta_v,'b-','LineWidth',2);
ax=gca;
ax.FontSize=20;
xlabel('\pi');
ylabel('\theta');
grid;
%f=gcf;
%exportgraphics(f,'fig_theta.png','Resolution',300);

figure(9),
plot(piv,w_v,'b-','LineWidth',2);
ax=gca;
ax.FontSize=20;
xlabel('\pi');
ylabel('Stationary wage rate $\bar{w}$','interpreter','latex');
grid;

figure(10),
plot(piv,pb_v(1,:),'b-',piv,pb_v(2,:),'k--','LineWidth',2);
ax=gca;
ax.FontSize=20;
xlabel('\pi');
ylabel('Lower threshold $\underline{p}_a$','interpreter','latex');
legend('$\underline{p}_a$ for $a=a_L$','$\underline{p}_a$ for $a=a_H$','interpreter','latex')
grid;

%% compute g'(0)
pi=0.03;
get_equilibrium;
dg0=0;
for i=1:2
    dg0 = dg0 + 0.5* (log(q(i)))^2 /vp(q(i),mu/pi)*mu/pi;
end
disp('(dg/ds)(0)');
disp(dg0); 